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ABSTRACT 

Context. Vortices in protoplanetary disks can capture solid particles and form planetary cores within shorter timescales than those 
involved in the standard core-accretion model. 

Aims. We investigate vortex generation in thin unmagnetized protoplanetary disks with an embedded giant planet with planet to star 
mass ratio 10"* and 10"'. 

Methods. Two-dimensional hydrodynamical simulations of a protoplanetary disk with a planet are performed using two different 
numerical methods. The results of the non-linear simulations are compared with a time-resolved modal analysis of the azimuthally 
averaged surface density profiles using linear perturbation theory. 

Results. Finite-difference methods implemented in polar coordinates generate vortices moving along the gap created by Neptune- 
mass to Jupiter-mass planets. The modal analysis shows that unstable modes are generated with growth rate of order 0.3flK for 
azimuthal numbers m = 4, 5, 6, where is the local Keplerian frequency. Shock-capturing Cartesian-grid codes do not generate very 
much vorticity around a giant planet in a standard protoplanetary disk. Modal calculations confirm that the obtained radial profiles of 
density are less susceptible to the growth of linear modes on timescales of several hundreds of orbital periods. Navier-Stokes viscosity 
of the order v = 10"^ (in units of a^Clp) is found to have a stabilizing effect and prevents the formation of vortices. This result holds 
at high resolution runs and using different types of boundary conditions. 

Conclusions. Giant protoplanets of Neptune-mass to Jupiter-mass can excite the Rossby wave instability and generate vortices in 
thin disks. The presence of vortices in protoplanetary disks has implications for planet formation, orbital migration, and angular 
momentum transport in disks. 

Key words. Planet and satellites: general - Accretion, accretion disks - fJydrodynamics - Instabilities - Methods: numerical 



1. Introduction 

Stability of rotationally supported gas disks is an area of active 
research, motivated among other reasons, by a need to under- 
stand the origin and stability of hydrodynamics turbulence un- 
derlying the so-called anomalous viscosity in accretion disks. 
The concept of a-turbulence in accretion disks was introduced 
more than three decades ago by Sha kura & Sunyaev| ( |1973| l to 
account for the angular momentum transfer and explain accre- 
tion onto the central object. The magnetorotational instability 
(MRI) has been proposed to explain the enhanced viscosity in 
hot and sufficiently ionized accretion disks with a Keplerian an- 
gular velocity profile threaded by a weak magnetic field (|Balbus| 

1998| l. 



& Hawley 1991 Balbus et al.| 1996 Balbus & Hawley 



However, in the context of cold protoplanetary disks, the ioniza- 
tion by cosmic rays and stellar radiation is limited to the surface 
layers of the disk while the so called "dead zone" in the vicinity 
of the central plane is expected to have low ionization ( Gammie 
|1996 V In some astrophysical systems such as cataclysmic vari- 
ables and outer regions of active galactic nuclei the coupling be- 
tween the magnetic field and the gas is also weak and MHD 
effects may be negligible. 
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The stability of diiTerentially rotating disks has been con- 
sidered analytically and numerically in the purely hydrodynam- 



ical case (Papaloizou & Pringle 



1984 1985 Goldreich et al. 



|1986[ |Papaloizou & Pringle| 1987 1) with applications to circum- 
stellar disks and galactic disks. A rotating isentropic torus with 
a gradient of specific angular momentum is found to be un- 
stable to low-order non-axisymmetric perturbations due to the 
Papaloizou-Pringle instability. Several mechanisms have been 
proposed that are able to sustain purely hydrodynamical turbu- 
lence and generate an anomalous a-viscosity in accretion disks 
(Li et al. 2000; Klahr & Bodenheimer 2003; Mukhopadhya^ 
iet al.,2005^ . |Dubrulle et al., ( ,2005 ) studied non-axisymmetric in- 
stabilities in stratified Keplerian disks using numerical and ana- 
lytical methods. A linear instability appears for Reynolds num- 
bers of order 10^ and perturbations with characteristic scales 
smaller than the vertical scale of the disk, assuming the angular 
velocity decreases with radius. These results suggest that despite 
the stabilizing eff'ect of the Coriolis force, a Keplerian flow may 
undergo a transition to turbulence. Nevertheless, some of those 
mechanisms may depend on boundary or edge effects. 

Rossby waves in thin Keplerian disks have been studied in 



the linear approximation ( [Lovelace et al. [1999[ [Li e t al."2000|l 
and with fully non-linear numerical simulations ([ Tagger_200l] l. 
The existence of unstable modes has been found to be associated 
with radial gradients of an entropy-modified version of vorten- 
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sity. Rossby waves in disks break up forming vortices in the non- 
linear limit dLi et al. 2001 1 in agreement with the predictions 
from linear theory. The dispersion relation of this Rossby Wave 
Instability (RWI) is analogous to the one for Rossby waves in 
planetary atmospheres. Reynolds stresses produced by the RWI 
can yield outward transport of angular momentum and accre- 
tion onto the central star. Varniere & Tagger ( 2006 1 have dis- 
cussed the generation of Rossby waves in the "dead zones" of 
protoplanetary disks where they may enhance the accretion rate 
of solids and favor planet formation. The angular momentum 
transport in disks around supermassive black holes at the center 
of galaxies can also be explained by the formation of Rossby 
vortices when there is a steep enough density gradient ( Colgate^ 
et al.|200? i. The angular momentum transfer due to vortices in 



galactic disks is found to be greater than in an a-viscosity disk. 
Recently, Tagger & Melia (2006 ) have described a magnetohy- 
drodynamics version of the RWI and applied it to the study of the 
quasiperiodic oscillations in Sgr A*. Rossby waves can also ap- 
pear in thin planetary atmospheres with solid rotation leading to 



the formation of vortices like Jupiter's "Great red spot" ( [Marcus 
[T988] l. 

Long-lived vortices are able to capture solid materials to 
form massive bodies and speed up the formation of planetary 



cores (|Barge & Sommeria 1995 Bracco et al. I999[ Klahr 



|& Bodenheimer| 2006| l. The stability of three-dimensional vor 



tices in a three dimensional stratified disk has been studied by 



Barranco & Marcus (2005 , 2006) using spectral anelastic hydro- 



dynamics simulations. They find that vortices are hydrodynam- 
ically stable for several orbits away from the mid-plane of the 
disk. The formation of vortices in the corotation region excited 
by a protoplanet has been studied numerically by Balmforth & 



Korycansky (20011 including the saturation of the corotation 



torque and the efl'ects of dissipation in the non-linear dynamics 
of the flow. 



Several numerical schemes studied by de Val-Borro et al. 
( [200 6 ) show vortex formation but the aim of the study did not 
include a detailed investigation of vortex generation. In this pa- 
per, we intend to look at this formation process in more detail. 
We also wish to check that codes that do not predict vortex gen- 
eration do not artificially damp unstable modes due to the nu- 
merical viscosity. Some simulations produce waves and vortices 
at the edge of the gap, which could in principle interact with 
the wake to cause semi-periodic disturbances propagating away 
along the shock. In simulations using the Piecewise Parabolic 
Method (hereafter PPM), low-m perturbations are observed at 
the edges of the gap (Ciecielag et al. 2000). Wave -like distur- 
bances with mode number ~ 5 are observed at the edge of the 
gap created by a Jupiter-size planet in the numerical results pre- 
sented by |Nelsorr & Benz (2003 ), Instabilities close to the planet 
or along the edges of the gap created by a giant planet, as well 
as the time variability of the flow near the Roche lobe may affect 
the speed and direction of planetary migration. 

In this paper, we study the effect of an annular gap cleared 
by a planet on the stability of a protoplanetary disk. We con- 
sider non-axisymmetric linear perturbations to the inviscid and 
compressible Euler equations. In Section 2, we present the semi 
analytical methods used to study the stability of disks. We de- 
scribe the numerical codes in Section 3. In Section 4 we present 
the results of the numerical simulations and the perturbative lin- 
ear analysis. We discuss the results in the context of protoplan- 
etary disks in Section 5. Finally, the numerical diffusivity in our 
numerical codes is calibrated in Appendix [A| 



2. Modal analysis 

We perform a modal analysis of analytical and numerically ob- 
tained density profiles, in order to see if there is agreement be- 
tween the vortex generation in the simulations and growing un- 
stable modes in the linear stability analysis. Linear perturbative 
analysis provides a valuable tool to study the stability of disks. 
The solution of the linearized Euler equations is treated along the 
lines of the work of [Lovelace et al.| ( fT999] ) and [Li et al.1 ( [2^0| ), 
where the stability can be evaluated solving a numerical eigen- 
value problem for a given profile. Alternatively, the growth of the 
initial perturbations can be determined by solving the equations 
as an initial value problem. 

We consider non-axisymmetric small perturbations sinu- 
soidally varying in azimuth to the inviscid Euler equations 



— + V ■ (£v) = 
ot 

d\ 1 - 

— -H (v ■ V)v = --VP - VO 
5f 2 



(1) 
(2) 



where £ is the surface mass density, P the vertically integrated 
pressure, v the velocity of the flow and <1> is the gravitational po- 
tential. The perturbed quantities have the form £ = Y.+SI.{r, (p, f), 
P - P + 5P{r, (p, t) and v = v -i- 6y(r, (p, t) with perturba- 
tions in the disk plane depending on the azimuthal angle as 
oc exp [i{m(p - (jDt)], where m is the azimuthal mode number and 
(jj - (ji)r + iy is the complex mode frequency. The initial velocity 
component in the radial direction is neglected and the angular 
velocity Q is obtained from the force balance between gravity, 
pressure gradients and centrifugal force 



1 



1 dP dO 
S dr dr 



(3) 



We use a barotropic equation of state, P oc with adiabatic 
index F = 5/3, to perform the modal analysis on the analytical 
disk profiles studied by [Lovelace et al. ( 1999| ). 

A locally isothermal equations of state is used to calculate 
unstable modes from the numerical simulations described in 
Section|3] The vertically integrated pressure is given by 



P = c% 



(4) 



where Cj is the sound speed. The temperature can be calculated 
from the ideal gas equation of state 



(5) 



where H is the gas constant and is the mean atomic weight. The 
initial unperturbed density is uniform in our simulations and the 
temperature has a fixed profile oc r ' for the isothermal calcula- 
tions, where r is the distance from the rotation axis. The sound 
speed profile of the disk is that of a standard slightly flaring solar 
nebula with constant aspect ratio 



c, = 0.05 



(6) 



where G is the gravitational constant and M» the mass of the 
central star. 

The linearized equations can be reduced to a second order 
differential equation for the enthalpy of the fluid, rj - 6P/'L, in 
the general case when the pressure is a function of both density 
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and temperature (|Loyelace et al.|[T999| |Li et al.|[2000l |de Val| 
Borro & Artymowicz|2004[ ). 



77" + B(r)j]' + C(r)j] = 



(7) 



where the coefficients B(r) and C(r) depend on the radial dis- 
tance. 



1 T' Q' 



1 - l; B(r) + Ak^^l|^u kyjAoj^ - 1 



and 

r(r) 



^2 _ _ cl/(L,Lp) 
m 



1 dZ2 
dr 



(8) 
(9) 

(10) 

(11) 
(12) 

(13) 
(14) 



Aa)(r) — LOr + iy — mD.{r) 

where k is the epicyclic frequency, / = Q.r^ is the angular mo- 
mentum per unit mass, Cj the adiabatic sound speed and Lj, 
Lp are the equilibrium length scales of variation of entropy and 
pressure changes, given by: 



1 d\n{PIYF) 
r 

1 dln(P) 
f dr 



(15) 
(16) 



Lj and Lp are calculated numerically from the averaged profiles 
obtained in the simulations using Equation]?] 

The growth of the unstable modes can form vortices or 
Rossby waves in the nonlinear regime ( |Li et al. 2001| l. The 
Rayleigh criterion states that the disk will be stable to axisym- 
metric perturbations when the specific angular momentum in- 
creases with radial distance. For a Keplerian disk the epicyclic 
frequency is always positive, - Q^, where the Keplerian an- 
gular frequency is given by Qk = (GMJry^, and therefore 
axisymmetric waves will be stable. However, when the pressure 
effects are taken into account it is possible to have axisymmetric 
instabilities for a sufficiently large pressure gradient according 
to the Solberg-H0iland criterion, 

K\r) + N^(r) > (17) 



where 



, 1 dP d / ^ \ 

^^^^^-fzd^d^^^y 



(18) 



is the square of the Brunt- Vaisala frequency in the radial direc- 
tion. 

The eigenproblem for the perturbed enthalpy is solved us- 
ing two semi-analytical methods. One way of finding the com- 
plex eigenfrequency is the shooting method, where the integra- 
tion proceeds from the disk boundaries to an intermediate fitting 
point, where continuity of the eigenfunction and its first deriva- 
tive is required. Our implementation uses a leapfrog method to 
integrate the equation from the boundaries to the fitting point. 
The values of the entalphy and its derivative at the starting points 
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Fig. 1. Real frequency of the most unstable modes for a gap 
opened by a Jupiter-mass planet after 10 periods as a function 
of the azimuthal mode number The solid line shows the mode 
frequency at the outer edge of the gap and the dashed line is 
the mode frequency at the inner edge divided by mOK('"o) where 
^K{ro) is the Keplerian angular frequency at the planet position. 
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Fig. 2. Growth rate for a gap opened by a Jupiter-mass planet 
after 10 orbits against the azimuthal mode number The solid 
line shows the growth rate at the outer edge of the gap and the 
dashed line is the growth rate at the inner edge. The growth rate 
of the instability peaks at mode numbers 5-6. 



are specified based on several prescriptions using outgoing spiral 
waves ( |Li et a l. 2000 ) and vanishing eigenfuntion at the bound- 
aries. We find that the obtained growth rates do not depend sensi- 
tively on the choice of boundary conditions. Several root finding 
algorithms in the complex plane can be used to find the unstable 
modes. The winding number theorem uses closed-path integrals 
(e.g. Ka rgl & Marston|1989 1 to find the number of roots inside 
the contour. 

The winding number theorem states that for a complex ana- 
lytic function /(w) defined inside a contour C 



2m(N - P) 



/V) 



(19) 
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where is the number of roots and P is the number of poles 
inside the contour C, and the integral is defined in counterclock- 
wise sense. We look for contours that comprise a single root and 
are not located close to solutions or singularities on the real axis. 
There are singularities in equation |7] at the corotation resonance 
Au) - and when a modified form of the Lindblad resonance 
for non-baro tropic flow is satisfied, - Act/ - c^/{LsLp) = 0. 
It is important to avoid branch points close to the contours since 
the winding number method would not give the right number of 
roots. In the case of analytical density profiles given by hyper- 
bolic trigonometric functions the branch points do not aff'ect the 
contours used to find solutions with positive growth rate. 
The solutions are then calculated by 



M 



!=1 



dw 



(20) 



where ^, are the zeros and the poles of the complex function 
/(w). The mode frequency can be determined in a contour with a 
single root. A multidimensional Newton-Raphson method (PressJ 
et al.||199"2 l is then employed to locate the roots with arbitrary 



accuracy. 

Another approach to solve equation [7] involves discretizing 
the equation on a finite grid and use appropriate boundary con- 
ditions to reduce the problem to finding numerically the roots of 



the determinant of a complex tridiagonal matrix (e.g., Laughlin 
et al.|1998 Li et al.|2000 i. We solved the determinant using the 



previous root finding algorithm and obtain the radial profile for 
the eigenfunction rjir) and the perturbed variables. 

We checked these two methods on the axisymmetric ana- 



lytical step jump profiles in surface density studied by |Li et al 
( |2000| l. We considered azimuthal mode numbers from m - 1 
to 10 and calculated the growth rates of the unstable modes and 
the corresponding eigenfunctions. For analytical density profiles 
with various shapes both our methods agree with the results of 



Li et al. (2000j|. We tested the dependence of the solution on 



the shape of the pressure profile and aspect ratio of the disk. The 
eigenfunction for density profiles with a locally isothermal equa- 
tion of state is obtained by solving the discretized equation]?] 

In Fig. [T] we show the real part of the mode frequency as a 
function of the azimuthal number for the averaged density pro- 
files of a NIRVANA simulation (see Section [3]l at x - 
256 X 768 resolution after 10 orbits. Fig. |2] shows the growth 
rate as a function of the mode number for the same averaged 
density profile, after 10 orbital periods when the threshold for 
the excitation of the instability is reached. The density slope in 
the simulations is described using two parameters. The depth of 
the gap is calculated with respect to the unperturbed density in 
the inner and outer disk. The length scale over which the density 
varies is estimated to be the diff'erence between the local maxi- 
mum and minimum at both inner and outer disks. An analytical 
jump function is fitted to the averaged profiles. 

The local maximum at the planet position in the averaged 
profiles is not expected to change the growth rates significantly. 
In addition, most of the gap is clean, apart from the material 
close to the planet position and at Lagrangian points L4, and L5. 
In the simulations with an initial gap the averaged profiles do not 
have local maxima inside the gap after 100 orbits. 

The larger growth rates in Fig. |2] correspond to the more un- 
stable modes of equation |7] that will dominate the solution. The 
inner edge of the gap opened by a planet, shown by the dashed 
line, has higher growth rate than that of the outer edge, repre- 
sented by the solid line, at a given time. This difference can be 
due to the fact that the inner boundary is closer to the location 
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Fig. 3. Radial profile of the eigenfunctions for the outer edge of 
the gap at f = 10 periods and mode number m - 5. From top to 
bottom the pressure perturbation and radial and azimuthal per- 
turbed velocity components are shown. The dotted and dashed 
lines are the real and imaginary part of the eigenfunctions. The 
amplitude is shown by the solid line which peaks at the position 
of the edge for the eigenfunction of the perturbed pressure. 



of the instability in the inner disk than it is the outer boundary 
to the corresponding instability outside the planet's radius. The 
real frequency of the modes correspond to a radial location just 
outside the corotating region where vortices are formed in the 
simulations after about 10 orbits. The number of vortices in the 
numerical results is consistent with the growth rates peaking at 
m — 4-6. 

In Fig. [3] the real and imaginary parts of the radial eigen- 
function of the perturbed variables at the outer edge of the gap 
are plotted for azimuthal mode m - 5. The eigenfunction corre- 
sponds to a time when the gap becomes deep enough to generate 
modes with positive growth rate. The middle and bottom panel 
show the eigenfunctions of the perturbations of the velocity com- 
ponents. The radial eigenfunctions at the inner edge of the gap 
after 10 orbits for azimuthal mode m = 5 are plotted in Fig.|4] 



3. Numerical codes 

We performed 2-dimensional hydrodynamical simulations using 
two independent grid-based codes, implemented in cylindrical 
and Cartesian coordinates. The simulations were run on a uni- 
form grid for 100 orbital periods Different boundary conditions 
were tested to avoid reflection of waves at the boundaries. 
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Fig. 4. Radial profile of the eigenfunctions for the inner edge of 
the gap at f = 10 periods and mode number m = 5. From top to 
bottom the pressure perturbation and radial and azimuthal per- 
turbed velocity components are shown. The dotted and dashed 
lines are the real and imaginary part of the eigenfunctions. 



3.1. Initial setup 

The computations were performed in the radial domain 0.4a < 
r < 2.5a where a is the planet semi-major axis. The disk is as- 
sumed to be geometrically thin and the fluid equations are solved 
using vertically-integrated variables: 



-f 

J-H 



p(r, 4>, z) dz, P{r, 



■, 0) = I pir, ( 

J-H 



,z)dz (21) 



where H is the vertical scale height. 

The planet's gravitational potential was given by the formula 

-1^ 



where e is the gravitational softening and s is the distance from 
the planet. The softening is defined as e = 0.6//p, with Hp the 
disc scale height at the planet location. This softening was in- 
troduced to reproduce the torque cut-ofi" due to 3-dimensional 
effects at distance from the planet larger than Hp. The soften- 
ing length is comparable in size to the characteristic size of the 
Roche lobe (along the planet's orbit) and thus it is likely to af- 
fect the flow dynamics in the corotation region but it is unlikely 
to strongly affect the dynamics in the regions where we observe 
vortex formation. Initially, the gas rotates with Keplerian angu- 
lar frequency around the central star The mass ratio had the 
values ju - Mp/M» - 10""' and 10""^, corresponding to Jupiter 
and Neptune masses when the stellar mass is one solar mass. 



Self-gravity of the disk was not included in the simulations. The 
planet was kept in a coplanar, circular orbit at a = 1 . We did not 
consider the accretion of disk material onto the protoplanet. 

We assume that the disk radiates efficiently the thermal en- 
ergy generated from tidal dissipation and viscous heating. In the 
absence of a radiation mechanism the gas would heat up and the 
disk could become geometrically thick. In our models, we use a 
locally isothermal equations of state, ( Lin & Papaloizou, 1985) 
with constant aspect ratio H/r = 0.05. 

The vortensity or potential vorticity is calculated in the coro- 
tating frame: 



f = (V X V 2np)/E 



(23) 



where Qp is the orbital frequency of the planet. 

We perform our simulations in the inviscid limit. In addi- 
tion, a few cases with physical viscosity v = 10"^ and 10"^ a^Q.^ 
are considered in order to estimate a threshold for the forma- 
tion of vortices. Artificial viscosity is not needed to smooth the 
shocks in our codes, which is important to study the formation 
of vortices. The codes are also able to resolve the large density 
contrast between the corotating region and the rest of the disk. 
Diffusion into the gap opened by the planet may result from nu- 
merical viscosity or shocks. The numerical viscosity in our codes 
is estimated in Appendix |A] for inviscid runs without tidal per- 
turbations. These tests confirm that numerical diffusion is not 
affecting our results. 

The initial density was uniform and the gravity of the planet 
was introduced gradually with the formula 



(TTt 

lOPf 



(24) 



We have used timescales between 5-10 orbital periods to intro- 
duce the gravity from the planet. Although the time when the 
instability appears depend on how the gravity is started, there 
is agreement between the modal analysis and numerical simula- 
tions for timescales of 5 and 10 orbits. The results presented in 
Section|4]use a switch-on time of the gravity of 5 orbits. 

In some of the calculations we use an initial gap profile de- 
rived under the WKB approximation (e.g., Lubow & D'Angelo 
|2006| ) 



I-oir) = exp 





3 











(25) 



where we use f/vp between 1.7 x 10^ and 1.8 x 10^ and Ap is 
the maximum of H and \r - a\. This gap profile allows us to 
check whether a rapid formation of the gap is necessary for the 
presence of instabilities. 

The Cartesian implementation of FLASH was run on a uni- 
(22) form non-rotating grid at resolution Wx x «y = 320 x 320, and 
Hx X Hy = 640 X 640. The computational domain was -2.6fl < 
X < 2.6a and -2.6ci <y< 2.6a. 

For computational convenience the unit of time used in the 
simulations was the orbital period at the planet location, Pp = 
27r[a^/G{M, + Mp)]'/^ = 2;r, where G(M, 4- Mp) = 1 and a = 1. 
Therefore, the angular frequency of the planet was Qp = 1 in our 
units. 

We used solid boundaries with wave damping zones next 
to the boundaries to avoid wave reflection that can create arti- 
ficial resonances. The damping regions were implemented in the 
NIRVANA simulations at 0.4fl < r < 0.5a and 2. la < r < 2.5a, 
by solving the following equation after each timestep: 

dx 

df " " 



R(r) 



(26) 
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where x represents the surface density and velocity components, 
T is the orbital period at the corresponding boundary, and R{r) 
is a parabolic function which is one at the domain boundary 
and zero at the interior boundary of the wave killing regions. 
NIRVANA simulations were also run using alternative boundary 
conditions to check that the gap profiles and vorticity in the coro- 
tating region do not depend strongly on boundary conditions. 

3.2. NIRVANA 

Our NIRVANA implementation is based on the original ver- 
sion of the code by Ziegler & Yorke| ( 1997[ l. The Navier-Stokes 
equations are solved using a directional operator splitting up- 
wind scheme which is second-order accurate in space and semi- 
second-order in time. 

NIRVANA uses a staggered grid where scalar quantities are 
stored at the center of grid cells and vectors are stored at the cell 
boundaries. The frame is centered on the center of mass of the 
star-planet system and rotates with the same angular frequency 
as the planet orbits the central star. The code was run with both 
non-reflective boundary conditions described above (see also 
Godon|1996 1 and wave-killing zones close to the boundaries ( de 
Val-Borro et al.|2006 i. These two types of boundary conditions 
provide consistent results since in both cases wave reflection in 
the boundaries is reduced. A Courant number of 0.5 was used in 
the simulations. 

Two grid resolutions were considered in which the number 
of cells in the radial and azimuthal directions were riy x - 
256 X 768 and 512 x 1536, respectively, with uniform spacing in 
both dimensions. Therefore, the cells around the planet position 
were approximately square. 

3.3. FLASH 



resolution close to the central star, an additional level of refine- 
ment was used in the inner disk in the Cartesian implementa- 
tion. Different timesteps were used for each level of refinement 
to speed up the simulation. 



The FLASH code ( |Fryxe U et aLpOOO)! is a fuUy parallel AMR 
implementati on of the PPM algorith m in its original Eulerian 
fomf] ( .Woodward & ColeUa.l984.,Coieria & Woodward! 1 984 . 



FLASH has been extensively tested in various compressible flow 



problems with astrophysical applications (Calder et al. 2002 
|Weirsetal.|2005] l. 

Our implementations of the FLASH code used both polar 
coordinates and the original Cartesian formulation of FLASH. 
The polar version of FLASH was run in the corotating coordi- 
nate system while the Cartesian implementation was run in the 
inertial frame. Our implementations were based on release 2.5 
of FLASH with customized modules for the equation of state 
and gravity forces, explicitly ensuring the conservative transport 
of angular momentum in the angular sweep. This is particularly 
important when large density gradients are present in the disk. 
The Coriolis forces were treated conservatively as described by 
|Kley| ( |1998| l. The isothermal Riemann solver was ported from 
the AMRA code ( |Plewa & MuUer|[200T] ). Courant numbers of 
0.7 and 0.8 were used in the simulations. 

The Cartesian grid was centered in the center of mass of 
the system and the orbits of the planet and the central star 
were integrated using a simple Runge-Kutta method. The grid 
cells were sized to give the same radial resolution in both im- 
plementations, although since the grid went to r = in the 
Cartesian simulations, the grid size had to be larger to achieve 
the same resolution. The Cartesian code was run at resolutions 
«x X ny = 320 X 320 and 640 x 640. To improve the angular 



The damping condition described in de Val-Borro et al. 
( |2006[ l was applied in the ring 2.1a < r < 2.5a close to the outer 
boundary but not in the inner disk for the Cartesian FLASH im- 
plementation. A free outflowing boundary was used at the outer 
boundaries x = ±2. 6a and y = ±2.6a. and there was free gas 
flow inside 0.4a. 

The cylindrical implementation used a frame centered on the 
star with resolutions rir x - 256 x 768 and 512 x 1536 and 
wave-damping zones close to the inner and outer boundaries. 
The Coriolis force, centrifugal force, and indirect terms due to 
the fact that the center of the frame is displaced from the center 
of mass of the system are included in the equation of motion. 

4. Results 

We carried out modal growth analysis on the density profiles 
of protoplanetary disks perturbed by an embedded giant planet 
ranging from a Neptune to a Jupiter mass. The time resolved 
linear analysis was performed on the azimuthally averaged den- 
sity profiles obtained from calculations, using a locally isother- 
mal equation of state. In the following description, we will show 
the results of inviscid runs at different resolutions with an em- 
bedded Jupiter-mass protoplanet using different boundary con- 
ditions. Results from runs with physical viscosity will also be 
presented. 

The real frequency and growth rate of the most unstable 
modes, as a function of time, are plotted in Fig.[5]for NIRVANA 
and FLASH polar simulations with resolution n,.x«0 = 256x768 
and wave-damping boundary conditions. The threshold for the 
appearance of the instability occurs after about 5 periods when 
the gap is sufficiently deep and the growth rate becomes pos- 
itive. The edge of the gap in the simulations becomes non- 
axisymmetric at this time and depressions in vortensity appear 
along the gap with m = 5-6 symmetry which coincide with the 
azimuthal number of the most unstable modes shown in Fig. |2] 
Small vortices are observed in the inner and outer disk shortly 
afterwards with the same angular distribution. 

In the bottom panel of Fig. |5] the growth rate from the inner 
disk edge, represented by crosses and circles for the different 
models, is larger than the growth rate at the outer disk, repre- 
sented by dots and stars. This difference may be artificial be- 
cause the inner boundary of the grid is closer to the planet loca- 
tion that the outer boundary is. The mode frequencies and growth 
rates for NIRVANA and polar FLASH agree within 10%. This 
is in agreement with the vortex sizes obtained from the Fourier 
analysis of the gravitational torques on the planet. In the end of 
the simulation, the growth rate of the instability is a fraction of 
the angular velocity at the edge of the gap. The disk becomes un- 
stable to axisymmetric perturbations at time ~ 40 orbits for the 



Jupiter simulations according to the Rayleigh criterion (Eq. 17 



FLASH is available at |http: //www.flash.uchicago.edu/ 



The RWI grows exponentially in agreement with the linear ana 
ysis during the first orbits ( Li et al.|20 0n and produces vortices 
that can be sustained by interaction with the spiral arms gener- 
ated by the planet. 

The linear analysis of Cartesian FLASH averaged profiles 
shows the presence of unstable modes in the inner and outer disk 
with growth rates of ~ 0.2Qk(/"o)- In some Cartesian FLASH 
runs there are mode solutions with large growth rates that appear 
at late times, when the gap is becoming deeper. At the same time 
there are indications that some mild instability is happening in 
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Fig. 5. Real frequency and growth rates of the unstable modes 
with azimuthal number m - 5, as a function of time, for 
NIRVANA and polar FLASH simulations with resolution rir x 
- 256 X 768. In the top panel, the crosses represent the mode 
frequencies at the outer edge of the gap and the dots are the mode 
frequencies at the inner edge for NIRVANA. The circles are the 
mode frequencies at the outer edge of the gap and the stars are 
the mode frequencies at the inner edge for polar FLASH. The 
frequencies are divided by mQK(ro), where Qk is the Keplerian 
angular frequency and ro is the radius at the edge of the gap. 
The crosses in the bottom panel are the growth rates at the outer 
gap edge and the dots are the growth rates at the inner edge for 
a NIRVANA calculation. The circles are the growth rates at the 
outer edge of the gap and the stars are the growth rates at the 
inner edge for a polar FLASH calculation. Growth rates are di- 
vided by the Keplerian frequency at the edge of the gap. 



the disk. However, these are not the fast-growing modes found 
in the polar grid calculations with azimuthal number m - 4-6. 
We do not find stationary solutions with positive growth rates 
that remain for a time of order of the growth timescale in the 
Cartesian grid models. This is consistent with the fact that no 
vortices appear in the edges of the gap in these simulations after 
tens of orbital periods, as shown in Fig.|9] 

Fig.|6]shows the density distribution for Jupiter inviscid sim- 
ulations in logarithmic scale using different boundary conditions 
for NIRVANA and FLASH models. The left panels show the 



Fig. 7. Growth rates of the unstable modes with mode number 
m = 5, as a function of time divided by the local Keplerian fre- 
quency. Dots represent a NIRVANA calculation with resolution 
tir X = 256 X 768 and stars are obtained from a FLASH cal- 
culation with resolution nrXnj, = 128 x 384. 
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Fig. 8. Azimuthally averaged surface density profiles for the 
Jupiter simulations after 100 orbital periods. The solid line is the 
NIRVANA simulation with resolution x - 256 x 768 and 
damping wave boundary conditions and the short dashed line is 
the NIRVANA simulation with resolution n,. xn^ - 256 x 768 
and outgoing-wave boundary conditions. The profiles are aver- 
aged over 5 orbits. 



density contours at f = 50 orbital periods and the right panels 
after f = 100 orbits. There are two vortices moving along the 
outer edge with different phase velocities in all the simulations 
at f = 50 periods. Those vortices will merge and form a large 
vortex at a later time producing strong oscillations on the torque 
exerted on the planet. 

In Fig. [T] we show the growth rates with azimuthal number 
m - 5 for the outer edge of the gap using an initial density with a 
gap given by Equation|25] The gap shape tends to a steady state 
towards the end of the simulation. The growth rates agree within 
5% (by the end of the simulation) with the results obtained using 
an initial uniform density. 
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Fig. 6. Surface density contours for cylindrical simulations in logarithmic scale. From top to bottom, NIRVANA simulation using 
a damping wave region as described in Section |3] (see also de Val-Borro et al!|2006^ , NIRVANA simulation using outgoing -wave 
boundary conditions defined by Godon|p996| l and polar FLASH. The left panels show the density after 50 orbital periods and the 
right panels show the density after 100 orbits with the same color scale. The resolution is x = 256 x 768 and the planet is 
located at xxy - (1 , 0) a. There are two elongated vortices forming next to the outer gap edge at 50 orbits that move with different 
velocities. The vortices merge in both simulations and one large vortex is observed in the outer disk at 100 orbits. 
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Table 1. The frequency in the corotating frame and magnitude 
of the maxima of the PDS of the gravitational torques from the 
inner disk are shown for the Jupiter simulations. The values are 
sorted by the magnitude of the PDS at the maximum. The magni- 
tude of the PDS is related with the amplitude of the oscillations 
at a given frequency. Polar codes have larger maxima and the 
frequencies correspond roughly to the Keplerian frequencies at 
the gap edges. 



Code 


Frequency 


Amplitude 


NIRVANA 


0.568 


2.81 X 10^" 


Flash polar 


0.531 


8.96 X 10-'* 


Flash Cart. 


1.611 


2.19 X 10-^ 



Table 2. The frequency in the corotating frame and magnitude 
of the maxima of the PDS of the gravitational torques from the 
outer disk are shown for the Jupiter simulations. 



Code 


Frequency 


Amplitude 


NIRVANA 


0.389 


1.11 X 10-^ 


Flash polar 


0.353 


4.82 X lO-*- 


Flash Cart 


0.997 


3.65 X 10-* 



Fig. [8] shows the profiles averaged over the azimuth for the 
NIRVANA simulations using various boundary conditions. The 
slope of the gap is steeper at the inner edge than at the outer 
edge for the cylindrical schemes. The density peak in the inner 
disk is about 50% of the peak at the outer gap edge. The over- 
all jump in density is larger at the inner edge despite the fact 
that the gap is slightly deeper just outside the planet's orbit. This 
produces a larger growth rate inside the corotation region (see 
Fig. [2]i. In the cylindrical FLASH simulation, the mass loss at 
the boundaries is greater but this does not affect the results of the 
modal calculation. The Cartesian FLASH results have a deeper 
gap and smaller density peaks at the edge of the gap. However, 
the growth rate of their most unstable modes is about 50% com- 
pared with that provided by the polar FLASH simulations. 

The size of the peaks in the power spectrum of the torques 
on a Jupiter-mass planet, from different regions in the disks, are 
shown in Tables [T] and |2] The torques are calculated excluding 
the material inside the Roche lobe of the planet, where the res- 
olution may not be good enough to resolve the circumplanetary 
disk. NIRVANA calculations have a larger peak, which corre- 
spond to the amplitude of the vortices moving along the edge of 
the gap. This is consistent with larger vortices being observed 
in the density distributions of NIRVANA calculations in Fig. [6] 
FLASH calculations have PDS peaks which are between one 
and two orders of magnitude smaller. The frequencies in the 
corotating frame are close to the local Keplerian angular fre- 
quency at the edges of the gap. The difference in the peak am- 
plitudes agrees with the results from the upwind and Godunov 
schemes studied by |de Val-Borro et al. ( j2006[ ). This correlation 
suggests that vortices can be formed by the non-linear evolution 
of Rossby waves in protoplanetary disks. 

The evolution of the vortensity in the NIRVANA simulation 
with resolution «i- x = 512 x 1536 is shown at several times 
in Fig. [To] The vortensity and Bernoulli constant are conserved 
for a barotropic inviscid fluid in the absence of discontinuities in 
the flow. In our case, the vortensity is roughly conserved along 
the streamlines in regions outside the Hill radius of the planet. 
Shock dissipation close to the protoplanet can lead to vortensity 
generation. As the gap is opened and strong trailing shocks are 
formed, the vortensity grows at the edge of the gap and along 



Table 3. The frequency and magnitude of the maxima of the PDS 
of the torques from the inner disk for the Neptune case are shown 
sorted by the magnitude of the PDS at the maximum. NIRVANA 
calculations have larger maxima with frequencies close to the 
Keplerian frequencies at the outer gap edge. 



Code 


Frequency 


Amplitude 


NIRVANA 


0.204 


4.67 X lO-'-" 


Flash polar 


0.221 


3.57 X 10-* 


Flash Cart. 


0.345 


3.99 X lO-** 



Table 4. The frequency and magnitude of the maxima of the 
PDS of the torques from the outer disk for the Neptune case are 
shown sorted by the magnitude of the PDS at the maximum. 



Code 


Frequency 


Amplitude 


NIRVANA 


0.389 


1.11 X 10-' 


Flash polar 


0.204 


5.27 X 10-** 


Flash Cart 


3.997 


5.411 X lO-'' 



the spiral arms. The vortensity at 10 orbital periods shows small 
cavities outside the peaks at the edge of the gap. These depres- 
sions break in 4-5 differentiated vortices when the growth rate of 
the RWI becomes positive. The vortensity peak at the outer gap 
edge and close to the outer spiral arm are corrugated while at the 
inner edge the vortensity is more stable. As explained before, 
this is probably an artificial effect due to the inner edge of the 
gap being closer to the inner boundary than the outer gap edge is 
to the outer boundary. The minima of vortensity rotating along 
the edge correspond to vortices observed in the density maps. In 
the bottom right panel of Fig. [TO] there is one single vortensity 
depression which is associated with a vortex located at azimuth 
a: n after 100 orbits. The vortensity inside the corotating region 
is considerably perturbed as the vortex moves along the gap. 

In Figs. 11 and 12 the azimuthally averaged vortensity in 
the inertial frame is shown. The initial vortensity profile for a 
disk with uniform density, has been subtracted. The 

Cartesian FLASH code has a large vortensity excess in the coro- 
tating region where the gap is more depleted than in the cylin- 
drical codes. The averaged vortensity in the outer disk is also 
greater in our Cartesian FLASH model. NIRVANA calculation 
shows vortensity peaks at the gap borders with a larger spike in 
the inner disk. 

The velocity fields plotted over the density contours in loga- 
rithmic scale at f = 100 orbits is shown in Fig. [13] The velocity 
vectors are calculated in the corotating frame of the local max- 
imum of pressure that coincides with the center of the vortex. 
The rotation in baroclinic sense is very clearly visible in the vor- 
tex close to the inner edge of the gap. The vortex in the outer 
disk interacts with the spiral wake created by the planet and is 
perturbed at this particular time. In Fig. 14 the streamlines are 
shown in the frame corotating with the vortex core. They show 
baroclinic rotation that is perturbed by interaction with the spiral 
wakes created by the planet. 



4.7. Dependence on physical viscosity 

In this Section we compare the unstable modes obtained from 
simulations including Navier-Stokes viscosity v = 10"^ and 10"^ 
(in code units, see Section[3]l. 

In Fig. 15 we show the growth rates in the outer disk of the 
unstable modes with m = 5 as a function of time for v - 10^^. 
The frequencies of unstable modes are calculated using perturba- 
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Fig. 9. Surface density distribution after 100 orbital periods for NIRVANA on the left hand side and FLASH on the right hand side 
using the same logarithmic color scale. Both models use the same wave damping condition in the outer disk between 2.1-2.5 a. 
while FLASH does not have a damping condition in the inner boundary. NIRVANA has density enhancements close to the gap 
opened by the protoplanet. FLASH has a smooth density distribution and a larger density peak at the planet position which is 
saturated in the image. 
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Fig. 11. Vortensity profiles averaged over azimuth at different 
times t = 10,20,50 and 100 orbital periods for the NIRVANA 
simulation at resolution x = 256 x 768. 



Fig. 12. Vortensity profiles averaged over azimuth at different 
times t - 10, 20, 50 and 100 orbital periods for the FLASH sim- 
ulation in Cartesian coordinates. 



tions on the inviscid Euler equations for simplicity. The growth 
rates at the outer edge of the gap agree within 20% with those 
obtained from inviscid simulations (see Fig.|5]l. However, notice 
that the linear analysis gives only an estimate of the growth rates 
in the viscous case. 



5. Discussion 

We have studied vortex formation in protoplanetary disks with 
an embedded giant planet, with mass ratios 10""* and lO""*, us- 
ing numerical simulations and linear perturbation analysis. The 
modal calculation is done following the strategy of |Lovelace| 



et al. ( |1999| l for a locally isothermal equation of state in a 
vertically-averaged disk. Vortices are formed in the cylindrical 
NIRVANA and FLASH 2-dimensional simulations in agreement 
with the linear analysis of non-axisymmetric perturbations. The 
growth rates calculated for NIRVANA and polar FLASH as a 
function of time agree within about 10%. The results of the lin- 
ear analysis are consistent with the absence of rapidly growing 
vortices near the edge of the gap in our Cartesian-grid PPM sim- 
ulations, which is thus not due to artificial numerical damping of 
unstable modes. This type of code does not produce the neces- 
sary steepness of the surface density profile and does not support 
growing non-axisymmetric perturbations. 
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Fig. 10. Vortensity in polar coordinates is shown at times f = 10, 20, 50 and 100 orbital periods from left to right and top to bottom. 
The simulation has resolution xn^ = 512 x 1536. 



Both our numerical schemes have very low numerical dif- 
fusion, which is estimated in Appendix |A] Runs with an ex- 
plicit Navier-Stokes viscosity were also performed. The unsta- 
ble modes in the outer disk calculated for v - 10"^ have growth 
rates ~ 0.25f2K at 100 orbits, which are 20% smaller than in the 
inviscid calculations. 

We speculate that the Cartesian-grid implementation may be 
more diffusive for this problem than polar geometry (in which 
the unperturbed Keplerian disk flows along the mesh structure) 
hence damping the growth of Rossby waves. The linear theory 
predicts that unstable modes will be present in the Cartesian sim- 
ulations after the gap is sufficiently deep but these modes have 
smaller growth rate than those obtained from NIRVANA and 
FLASH polar simulations. The growth rate for the inner edge 
of the gap is larger than it is at the outer gap edge. This may 



be artificial because the inner boundary is closer to the planet 
position that the outer boundary is to the outer gap edge. 

We observe a correlation between the growth rate of the un- 
stable modes in the linear analysis and the size of the peaks in 
the power spectrum of the gravitational torque on the planet by 
the disk. This correlation suggests that vortices in protoplane- 
tary disks can form close to the gap, produced by an embed- 
ded giant planet, from the collapse of Rossby waves. Vortices 
may grow and be sustained for long timescales by interaction 
with the planetary wake ( [Koller et al.|2003[p~et al.|2005 ). The 
two-dimensional approximation for the disk flow is anticipated 
to give qualitatively correct results although a three-dimensional 
analysis is needed to understand heat dissipation in the verti- 
cal direction and refraction effects in radially propagating waves 
( |Lin et al.|1990| l. An important restriction of our simulations is 
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Fig. 13. Velocity vectors in the corotating frame after 100 periods for the NIRVANA simulation at resolution x = 512 x 1536. 
The left panel shows the vortex at the inner gap edge and the right panel shows the vortex at the outer gap edge. 




that the planet is kept on a fixed circular orbit. It would be of 
interest to study how the vortices rotating along the edge of the 
gap affect the migration rate of a freely moving protoplanet em- 
bedded in a 3-dimensional disk. 

In summary, the linear analysis confirms that Rossby waves 
are formed in a thin protoplanetary disk with a giant planet with 
mass ratio between fi = lO^'^-lO"-', within tens of orbital pe- 
riods. The unstable modes with larger growth rates generate a 
non-axisymmetric perturbation, with mode number m - 4-6, 
which breaks into vortices in the nonlinear regime producing 
a non-axisymmetric density distribution. At the time when the 
growth rate becomes positive, small depressions in vortensity 
appear along the gap. These results do not depend on resolu- 
tion or boundary conditions. Simulations with an initial gap also 



generate vortices close to the edge of the gap. The growth rates 
estimated from the linear theory agree with the growth rates at 
late times in simulations with an initially flat density distribu- 
tion. A protoplanetary disk with a giant planet becomes popu- 
lated with vortices and spiral shocks that can efficiently trans- 
port angular momentum. This effect can be important in disks 
that are not sufficiently ionized to sustain turbulence via the MRI 
instability. We conclude that vorticity generation in protoplane- 
tary disks with an embedded giant planet is a robust mechanism 
that can lead to planet formation and radial transfer of angular 
momentum. 
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Fig. 15. Dependence on time of growth rates of modes with m - 
5 for NIRVANA simulations with resolution n^xn^ - 256 x 768 
(dots) and FLASH with resolution n^xn^ - 128 x 384 (stai's). 
The Navier-Stokes viscosity is v = 10"^. 
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Appendix A: Calibration of numerical viscosity 

The calculations presented in this paper are based on inviscid 
Euler equations. However, even though neither physical nor ar- 
tificial viscosity terms enter these equations, each numerical 
scheme has some intrinsic diffusivity that can be interpreted as a 
numerical viscosity, Vn. 

To calibrate the numerical viscosity in each of the hy- 
drodynamics codes, we use the numerical setup described in 
Section 3.1 with mass ratio q - Q. Tests are performed at 



two different resolutions to check for consistency of the results. 
Numerical viscosity is bound to depend on the flow properties. 
The values reported in this Appendix apply to disks in quiescent 
conditions and may therefore represent lower limits for numeri- 
cal diffusion in models with uniform initial density and fast gap 
formation. 

In the case of NIRVANA, a time-averaged measure of the 
numerical diffusivity is obtained over a 50 orbit period, as a 
function of the radial position, by analysing the trajectories of 
500 tracer (massless) particles released in the disk. The equa- 
tion of motion of each particle is integrated every hydrodynam- 
ics timestep by interpolating the velocity field at the particle's 
location and advancing in time its position by means of a second- 
order Runge-Kutta method. The spatial interpolation of the ve- 
locity field is also second-order accurate. Hence, trajectories are 
formally second-order accurate in both space and time. 

We assume that the viscosity is constant in a radial interval 
Ar - 0.1, which contains about 25 equally-spaced tracer parti- 
cles and is orders of magnitude larger than their diffusion length 
scale. We measure the averaged length travelled by the particles 
in each radial interval, over about 50 orbits (at r = 1), and es- 
timate the amount of numerical viscosity under the assumption 
that the particle drift velocity is 



3Vn 

2r ■ 



(A.l) 
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Experiments executed at resolutions n, x - 256 x 768 and 
512 X 1536 provide very similar results. 

The largest numerical diffusion is observed close to the i nner 



grid boundary, where Vn ~ 10 ^ (in code units, see Section 3.1 
at r Si 0.55 and ~ 10"^ atr^ 0.65 . In the radial domain between 
r =s 0.75 and r ~ 1.75, v„ lies between ~ 10"'° and ~ 10"''. In 
the outer part of the simulated disk, the numerical viscosity is 
comprised between ~ 10"^ and ~ 10"^. 

The numerical viscosity in FLASH is calibrated using a 2- 
dimensional local patch of a Keplerian disk with a massless sink 
hole in the center of the domain The sink hole has radius equal 
to the Roche radius 



1/3 

(A.2) 



for a protoplanet with mass ratio q - 10^^. The dimensions of 
the shearing box are 207?r x 20^r. We use local Cartesian coor- 
dinates in the Hill approximation corotating with the sink hole. 
The horizontal axis corresponds to the radial direction and the 
vertical axis to the direction of motion of the flow. Initially, the 
unperturbed surface density is uniform and the velocity of the 
matter flowing into the computational domain has only a verti- 
cal component given by the linearized Keplerian velocity 

3 

v.v = -i^^kx (A.3) 

We implement periodic boundaries in the vertical direction. Our 
units in the simulation are the Roche radius, 7?r, the initial sur- 
face density, Eo, and the Keplerian angular frequency. Ok, at the 
center of the shearing box. A cavity is produced at the orbital 
radius of the sink hole while the gas in the vicinity of the hole is 
pulled into it on the viscous timescale. 

The equation of viscous diffusion for a thin accretion disk 
can be obtained in the asymptotic limit assuming constant nu- 
merical diffusivity (Pringle 1981 ; Bry den et al.|1999 1. Using the 



OUti 



boundary conditions 2 = at x = 0, and dl,/dx - at x = 
where Xout - 10/?r is the distance from the center of the shearing 
box to the outer radial boundary, the surface density distribution 
for positive x is given by 

The kinematic viscosity can be evaluated from the ratio of sur- 
face densities at different times 



44 



out 1 
= T— 7- — In 



(A.5) 



We estimate values of Vn between 10"^ and 10"^ using the sur- 
face density profiles averaged in the vertical direction at times 
close to 50 orbits and grid resolutions of n^xtiy - 320 x 320 and 
640 X 640. 



